module discovery

The ModID function above is modified to accept a batch variable and to correct with combat. Combat was choosen as opposed to more single cell specific techniques (rPCA, Harmony, FastMNN, etc.) for a few reasons. 1. Low cell numbers - we don’t have enough different cells for our conditions for graph based methods to learn things efficiently. 2. Need corrected gene values - the above correct in a reduced PCA like space, not the gene expression values themselves. CCA is an alternative here but the concern from point 1 remains. 3. Reduce chance of overcorrection. Combat’s weakness in most applications is that it corrects only weekly and this is not sufficient to overcome tech or batch effects when applying unsupervised methods across those variables. It does not (at least in the studies I’ve seen) force disparate cell types together in the way that CCA can, which otherwise makes the DE testing have elevated false discovery.

results<-subset(results, idents = "Proliferating")
## Loading required package: Signac
results$batch<-case_when(grepl("4|5|6",results$orig.ident)~1,grepl("7|8|9",results$orig.ident)~2,grepl("10|11|12",results$orig.ident)~3 )
mods<-ModID(results, combat=TRUE, batch=results$batch, nPCS = 8, softpower = 16)
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at -2.4183
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.30103
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 9.8548e-16
## PC_ 1 
## Positive:  ANK3, TSHZ2, ZBTB20, FAAH2, MAML2, PRKCA, LDLRAD4, PCAT1, ABLIM1, TTN 
##     MGAT4A, C22orf34, UST, AL136456.1, PFKFB3, TAFA1, SESN3, MPP7, LINC01422, AL353660.1 
## Negative:  TUBA1B, TUBB, HIST1H4C, GAPDH, RRM2, STMN1, DIAPH3, ACTB, PFN1, HIST1H3B 
##     POLQ, TYMS, NCAPG, TK1, HIST1H2AH, TMPO, PCLAF, HIST1H1B, CFL1, HMGN2 
## PC_ 2 
## Positive:  POLQ, NCAPG, DIAPH3, SPC25, CIT, KIF18B, KIF15, RRM2, ASPM, ATAD5 
##     NCAPG2, CDCA2, NUSAP1, KIF11, KIF14, TOP2A, KNL1, BUB1B, KIF4A, BRIP1 
## Negative:  RPL28, RPL18A, RPL41, RPS15, RPL27A, RPL13A, RPS28, RPS19, RPL11, RPS12 
##     RPS23, RPS29, RPS18, RPS27, RPL10, RPL18, RPL37, RPS16, B2M, RPLP1 
## PC_ 3 
## Positive:  ABLIM1, ZBTB20, MGAT5, MAP3K1, RORA, WDFY2, SSBP2, MAP3K4, RFX3, MCF2L2 
##     TSHZ2, ARL15, AL136456.1, DUSP16, AKT3, ATP13A3, PHF21A, UST, MAST4, PRKCA 
## Negative:  MTRNR2L8, NKG7, AGAP1, PECAM1, CCL5, GZMA, MKI67, CD8A, CENPF, HMGB2 
##     KLRK1, CD8B, AOAH, ZEB2, KPNA2, NAALADL2, TOP2A, HIST1H1E, HSP90B1, INCENP 
## PC_ 4 
## Positive:  PTPN3, MIR4422HG, PECAM1, DTHD1, LYST, SESTD1, GZMK, NCOA1, MCTP2, CNR2 
##     HAVCR2, AL031599.1, UBE2E2, AGAP1, AFF3, KLRK1, ANTXR2, TRPS1, SAMD3, RPTOR 
## Negative:  LTB, SORL1, MAP3K1, ANK3, TTC39C, KIF14, TSHZ2, CDCA2, KIF4A, CKAP2L 
##     KIF2C, TNFAIP3, CDCA8, JUNB, UBE2C, ESPL1, RPS19, EPHA4, FMN1, SGO2 
## PC_ 5 
## Positive:  DTL, MSH6, CLSPN, TRERF1, EXO1, CDC45, FANCA, MCM4, CDC6, CCNE2 
##     CTSW, HELLS, CCL5, ZEB2, POLE2, MCOLN2, C1orf21, MTHFD1L, FAM111B, THEMIS 
## Negative:  SESN3, SESTD1, MIR4422HG, KIF14, PTPN3, LDLRAD4, PLK1, KIF20A, KIF23, AFF3 
##     KIF4A, HMMR, BACH2, AURKB, TXK, DTHD1, CCNB2, APOLD1, AURKA, CDCA8 
## Found 7 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.
## Found3batches
## Adjusting for0covariate(s) or covariate level(s)
## Standardizing Data across genes
## Fitting L/S model and finding priors
## Finding parametric adjustments
## Adjusting the Data
##    Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1      1   0.0541  2.49         0.3810 157.000  157.0000 176.00
## 2      2   0.1100 -1.56         0.3420  85.300   83.7000 109.00
## 3      3   0.4000 -2.11         0.4900  47.400   45.0000  70.90
## 4      4   0.6220 -2.05         0.6440  27.100   24.4000  48.20
## 5      5   0.6670 -2.04         0.6290  16.000   13.3000  34.10
## 6      6   0.6800 -1.89         0.6630   9.730    7.3000  24.90
## 7      7   0.6400 -1.79         0.6040   6.150    4.1000  18.70
## 8      8   0.2200 -3.92         0.1170   4.020    2.2900  14.40
## 9      9   0.2210 -3.47         0.1360   2.710    1.3100  11.20
## 10    10   0.2240 -3.12         0.1440   1.890    0.8090   8.86
## 11    12   0.7490 -1.37         0.7940   0.983    0.3050   5.72
## 12    14   0.7670 -1.24         0.8080   0.554    0.1240   3.81
## 13    16   0.7710 -1.22         0.8000   0.330    0.0490   2.59
## 14    18   0.1990 -2.04         0.0217   0.204    0.0222   1.80
## 15    20   0.2000 -1.94         0.0272   0.130    0.0091   1.27

## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.

## dynamicMods
##   0   1   2   3   4   5   6 
## 212  21  15  14  13  12  11

##  mergeCloseModules: Merging modules whose distance is less than 0.5
##    multiSetMEs: Calculating module MEs.
##      Working on set 1 ...
##      moduleEigengenes: Calculating 7 module eigengenes in given set.
##    multiSetMEs: Calculating module MEs.
##      Working on set 1 ...
##      moduleEigengenes: Calculating 5 module eigengenes in given set.
##    Calculating new MEs...
##    multiSetMEs: Calculating module MEs.
##      Working on set 1 ...
##      moduleEigengenes: Calculating 5 module eigengenes in given set.

## $yellow
##  [1] "PCAT1"      "UST"        "MPP7"       "EPHA4"      "SSBP2"     
##  [6] "AOAH"       "SESTD1"     "JUNB"       "AP000787.1" "NDUFA6"    
## [11] "H2AFZ"      "PKM"        "DISC1"     
## 
## $red
##  [1] "TUBA1B"   "TUBB"     "DIAPH3"   "POLQ"     "NCAPG"    "PCLAF"   
##  [7] "HIST1H1B" "PCNA"     "ASPM"     "CIT"      "KIF15"    "ATAD5"   
## [13] "NCAPG2"   "KIF11"    "TOP2A"    "KNL1"     "BRIP1"    "CIP2A"   
## [19] "BRCA2"    "ECT2"     "FANCI"    "CENPF"    "DTL"      "CLSPN"   
## [25] "FANCA"    "HELLS"    "HNRNPAB"  "CENPP"    "C12orf75" "MCM5"    
## [31] "MCM7"     "KNTC1"   
## 
## $brown
##  [1] "STMN1"   "DUSP16"  "STAM"    "CDHR3"   "CCL5"    "HSP90B1" "MCTP2"  
##  [8] "ANTXR2"  "APOLD1"  "SCMH1"   "AGPAT4"  "ANXA2"   "SYTL2"   "TGFBR3" 
## 
## $blue
##  [1] "RPL28"  "RPL18A" "RPL41"  "RPS15"  "RPL27A" "RPL13A" "RPS28"  "RPS19" 
##  [9] "RPL11"  "RPS12"  "RPS23"  "RPS29"  "RPS18"  "RPS27"  "RPL10"  "RPL18" 
## [17] "RPL37"  "RPS16"  "B2M"    "RPLP1"  "RPS2"   "RPS27A" "RPL37A" "RPS15A"
## [25] "RPLP2"  "TMSB10" "RPS24" 
## 
## [1] "0/10000 permutations failed the Mann-Whitney test."
## [1] "0/10000 permutations failed the Mann-Whitney test."
## [1] "13/10000 permutations failed the Mann-Whitney test."
## [1] "0/10000 permutations failed the Mann-Whitney test."
## $yellow
##  [1] "PCAT1"      "UST"        "MPP7"       "EPHA4"      "SSBP2"     
##  [6] "AOAH"       "SESTD1"     "JUNB"       "AP000787.1" "NDUFA6"    
## [11] "H2AFZ"      "PKM"        "DISC1"     
## 
## $red
##  [1] "TUBA1B"   "TUBB"     "DIAPH3"   "POLQ"     "NCAPG"    "PCLAF"   
##  [7] "HIST1H1B" "PCNA"     "ASPM"     "CIT"      "KIF15"    "ATAD5"   
## [13] "NCAPG2"   "KIF11"    "TOP2A"    "KNL1"     "BRIP1"    "CIP2A"   
## [19] "BRCA2"    "ECT2"     "FANCI"    "CENPF"    "DTL"      "CLSPN"   
## [25] "FANCA"    "HELLS"    "HNRNPAB"  "CENPP"    "C12orf75" "MCM5"    
## [31] "MCM7"     "KNTC1"   
## 
## $brown
##  [1] "STMN1"   "DUSP16"  "STAM"    "CDHR3"   "CCL5"    "HSP90B1" "MCTP2"  
##  [8] "ANTXR2"  "APOLD1"  "SCMH1"   "AGPAT4"  "ANXA2"   "SYTL2"   "TGFBR3" 
## 
## $blue
##  [1] "RPL28"  "RPL18A" "RPL41"  "RPS15"  "RPL27A" "RPL13A" "RPS28"  "RPS19" 
##  [9] "RPL11"  "RPS12"  "RPS23"  "RPS29"  "RPS18"  "RPS27"  "RPL10"  "RPL18" 
## [17] "RPL37"  "RPS16"  "B2M"    "RPLP1"  "RPS2"   "RPS27A" "RPL37A" "RPS15A"
## [25] "RPLP2"  "TMSB10" "RPS24"
saveRDS(mods, "~/gibbs/DOGMAMORPH/Ranalysis/modules/20230601ProliferatingModules.rds")

for (i in names(mods)){
  write.table(mods[[i]], file=paste0("~/gibbs/DOGMAMORPH/Ranalysis/modules/proliferation/combat_",i,".txt"), quote = FALSE, row.names = FALSE, col.names = FALSE, sep = "")
}

annotating modules

After correction, we find 4 modules, of those four none appear paritcularly different across treatment, batch, or timepoint. Individual to individual variation is present, but its ultimately not aligned with an interesting variable.

Further iterations are below either doing the batch effect correction before feature identification or without merging modules based on similarity. In all, it looks like the only two consistent signal (at a module level) in proliferating cells is that: 1. Proliferation 2. Ribosomal/Housekeeping

results<-AddModuleScore(results, mods, name = names(mods))
modnames<-paste0(names(mods), seq(1, length(names(mods))))
FeaturePlot(results, features = modnames, min.cutoff = 'q5', max.cutoff = 'q95')

VlnPlot(results, features = modnames, split.by = "Participant")
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##       
## This message will be shown once per session.

VlnPlot(results, features = modnames, split.by = "batch")

VlnPlot(results, features = modnames, split.by = "Treatment")

VlnPlot(results, features = modnames, split.by = "Timepoint")

#borrowing some rushmore colors from the wes anderson color pack
BottleRocket2 = c("#FAD510", "#CB2314", "#273046", "#354823", "#1E1E1E")


#defining a quick function to read in and plot Enrichr 

PlotEnrichment<-function(file, topn=10, returnplot=TRUE, title="Significant Pathways"){
  data<-read.table(file, sep = "\t", header=1)
  data$value<- -log10(data$Adjusted.P.value)
  #order by most signfiicant 
  data$Term<-factor(data$Term, levels = rev(data$Term))
  data<-data[order(data$value, decreasing = TRUE),]
  #subset to top n 
  data<-data[1:topn,]
  plot<-ggplot(data, aes(x=value, y=Term))+geom_bar(stat="identity", fill=BottleRocket2[2])+xlab("-log10(p_adj)")+ylab("Pathway")+ggtitle(title)+theme_classic()
  if(returnplot){return(plot)}
  return(data)
}

todo<-list.files("~/gibbs/DOGMAMORPH/Ranalysis/modules/proliferation/enrichr_combat/")

for (i in todo){
  print(PlotEnrichment(paste0("~/gibbs/DOGMAMORPH/Ranalysis/modules/proliferation/enrichr_combat/",i), topn = 20, title=i))
}

running combat before running module ID

Rational here is that if we correct for batch before doing variable features and PCA we may be able to discover more true signals. In reality, we end up very consistently recovering proliferation and ribosomal modules only despite different npcs and softpowers.

ModID<-function(seuratobj, cluster=NA, prefix="modules", nPCS="", softpower="", combat=FALSE, batch=NA){
    #subset to a cluster or clusters of interest if desired 
    if (sum(!is.na(cluster))>1){test_clus<-subset(seuratobj, idents=cluster)}else{test_clus<-seuratobj}
    
    #perform feature selection, scaling, and PCA 
    if(combat==TRUE){
        test_clus_2<-test_clus@assays$RNA@data
        test_clus_2<-ComBat(test_clus_2,batch)
        test_clus@assays$RNA@data <- test_clus_2
    }
    test_clus<-FindVariableFeatures(test_clus, verbose=FALSE)
    test_clus<-ScaleData(test_clus, verbose=FALSE)
    test_clus<-RunPCA(test_clus, verbose=FALSE)
    
    #prompot for number of PCs needed 
    if(nPCS==""){print(ElbowPlot(test_clus))}
    while(nPCS==""){nPCS<-readline("how many NPCs?")}
    nPCS<-as.integer(nPCS)
    #expanding the genes by projecting the dims to more PCs 
    test_clus<-ProjectDim(test_clus)
    
    #gere we'll grab the top 25 features (positive and negative) associated with each dimension 
    genes<-c()
    for (i in 1:nPCS){genelist<-TopFeatures(test_clus, dim = i, nfeatures = 50,balanced = T, projected = T )
    for (i in 1:length(genelist)){genes<-c(genes, genelist[[i]])}}
    #remove duplicates 
    genes<-unique(genes)
    #grab the matrix from the seurat object 
    test_clus<-as.matrix(test_clus@assays$RNA@data[genes,])
    
    #prompt for softpower 
    FindPower(datExpr=t(test_clus))
    while(softpower==""){softpower<-readline("what softpower?")}
    softpower<-as.integer(softpower)
    #run the remaining functions to actually cut the tree
    test_clus_tom<-ClusterTOM(datExpr = t(test_clus), softPower = softpower)
    test_clus_mods<-CutTOMTree(datExpr = t(test_clus), geneTree = test_clus_tom$geneTree, dissTOM = test_clus_tom$dissTOM, minModuleSize = 10 )
    test_clus_merge_mods<-MergeSimilarModules(datExpr=t(test_clus), dynamicColors = test_clus_mods$dyanmicColors, geneTree = test_clus_tom$geneTree, MEDissThres = 0.5)
    print(test_clus_merge_mods$merged_modules)                                          
    test_clus_merge_mods.isSig = sapply(test_clus_merge_mods$merged_modules, function(module){
        TestModuleSignificance(mod = module, dissTOM = test_clus_tom$dissTOM, expr.data = test_clus,
                               n_perm = 10000, pval = 0.05, n.bin = 10)
    })
    test_clus_merge_mods.isSig = test_clus_merge_mods$merged_modules[test_clus_merge_mods.isSig]
    print(test_clus_merge_mods.isSig)
    names(test_clus_merge_mods.isSig)<-paste(prefix,names(test_clus_merge_mods.isSig ), sep="_")
    return(test_clus_merge_mods.isSig)
}

#tried npcs/softpower 
#8/18 - Ribosomal and a proliferating 
#7/18 - Ribosomal and a proliferating 
#5/18 - Ribosomal and a proliferating 
#8/10 - Ribosomal and a proliferating 
#8/5 - Ribosomal and a proliferating 
#5/5 - Ribosomal and a proliferating

mods<-ModID(results, combat = TRUE, batch = results$batch, nPCS = 8, softpower = 18)
## Found 24380 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.
## Found3batches
## Adjusting for0covariate(s) or covariate level(s)
## Standardizing Data across genes
## Fitting L/S model and finding priors
## Finding parametric adjustments
## Adjusting the Data
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at -2.4183
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.30103
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 9.8548e-16
## PC_ 1 
## Positive:  ANK3, ZBTB20, PCAT1, FAAH2, TSHZ2, MAML2, ABLIM1, PRKCA, SESN3, TTN 
##     LDLRAD4, UST, MPP7, FMN1, AL353660.1, ACYP2, TAFA1, MGAT4A, AL136456.1, EPHA4 
## Negative:  TUBA1B, TUBB, HIST1H4C, RRM2, STMN1, DIAPH3, HIST1H3B, GAPDH, POLQ, ACTB 
##     HIST1H2AH, NCAPG, TK1, HIST1H1B, PFN1, TMPO, NUSAP1, PCLAF, SHCBP1, ASPM 
## PC_ 2 
## Positive:  POLQ, NCAPG, SPC25, DIAPH3, CIT, KIF15, KIF18B, RRM2, ASPM, KIF11 
##     CDCA2, NCAPG2, ATAD5, NUSAP1, KIF14, KIF4A, TOP2A, BUB1B, KNL1, BRIP1 
## Negative:  RPL28, RPL18A, RPL41, RPL13A, RPS15, RPL27A, RPS19, RPL11, RPS12, RPS28 
##     RPS23, RPS29, RPS18, RPL10, RPL37, RPS27, RPS2, RPLP1, B2M, RPS16 
## PC_ 3 
## Positive:  ABLIM1, MGAT5, RORA, MAP3K1, ZBTB20, SSBP2, WDFY2, MAP3K4, MCF2L2, ATP13A3 
##     AL136456.1, RFX3, ARL15, TSHZ2, PHF21A, STAM, DUSP16, AKT3, MAST4, UST 
## Negative:  MTRNR2L8, NKG7, AGAP1, CCL5, PECAM1, CD8A, MKI67, CENPF, GZMA, KLRK1 
##     HMGB2, TOP2A, CD8B, INCENP, KPNA2, AOAH, HIST1H1E, CENPE, CDCA2, DLGAP5 
## PC_ 4 
## Positive:  PECAM1, PTPN3, DTHD1, MIR4422HG, LYST, MCTP2, NCOA1, GZMK, KLRK1, CNR2 
##     SESTD1, NKG7, AGAP1, AL031599.1, HAVCR2, KCNQ5, UBE2E2, SAMD3, GZMA, AFF3 
## Negative:  LTB, MAP3K1, SORL1, KIF14, KIF4A, CDCA2, CDCA8, UBE2C, KIF23, TSHZ2 
##     ANK3, CKAP2L, KIF2C, PLK1, KIF20A, CCNB2, ESPL1, CCNF, GTSE1, HMMR 
## PC_ 5 
## Positive:  MIR4422HG, PTPN3, SESTD1, DTHD1, AFF3, GZMK, SESN3, AL031599.1, LDLRAD4, PECAM1 
##     UBE2E2, ANTXR2, BACH2, LYST, HAVCR2, KIF20A, PLK1, DSCAML1, PACSIN1, KIF14 
## Negative:  EXO1, DTL, FANCA, CLSPN, TRERF1, CCNE2, CDC6, MSH6, CDC45, CTSW 
##     MCM4, POLE2, FAM111B, HELLS, MCOLN2, CENPP, C1orf21, CCL5, MPP6, MCM10 
##    Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1      1   0.0259  1.50          0.270 155.000  1.54e+02 173.00
## 2      2   0.1620 -1.85          0.310  83.900  8.23e+01 107.00
## 3      3   0.4210 -2.20          0.444  46.700  4.41e+01  70.00
## 4      4   0.6170 -2.07          0.595  26.800  2.38e+01  47.80
## 5      5   0.6470 -1.92          0.591  15.900  1.29e+01  34.00
## 6      6   0.6510 -1.81          0.615   9.800  7.19e+00  25.00
## 7      7   0.6210 -1.70          0.582   6.250  3.96e+00  18.90
## 8      8   0.6410 -1.54          0.599   4.140  2.28e+00  14.60
## 9      9   0.6930 -1.44          0.669   2.830  1.34e+00  11.50
## 10    10   0.6790 -1.33          0.653   1.990  8.09e-01   9.10
## 11    12   0.2150 -2.62          0.131   1.070  3.02e-01   5.91
## 12    14   0.2100 -2.26          0.143   0.612  1.22e-01   3.96
## 13    16   0.2080 -2.13          0.146   0.370  5.02e-02   2.71
## 14    18   0.8330 -1.12          0.877   0.231  2.27e-02   1.90
## 15    20   0.8060 -1.11          0.840   0.148  9.02e-03   1.35

## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.

## dynamicMods
##   0   1   2   3   4 
## 195  35  25  25  13

##  mergeCloseModules: Merging modules whose distance is less than 0.5
##    multiSetMEs: Calculating module MEs.
##      Working on set 1 ...
##      moduleEigengenes: Calculating 5 module eigengenes in given set.
##    multiSetMEs: Calculating module MEs.
##      Working on set 1 ...
##      moduleEigengenes: Calculating 4 module eigengenes in given set.
##    Calculating new MEs...
##    multiSetMEs: Calculating module MEs.
##      Working on set 1 ...
##      moduleEigengenes: Calculating 4 module eigengenes in given set.

## $brown
##  [1] "PCAT1"    "FAAH2"    "UST"      "ACYP2"    "EPHA4"    "CMTM8"   
##  [7] "ZNF331"   "HIST1H4C" "STAM"     "DUSP16"   "CDHR3"    "CCL5"    
## [13] "MCTP2"    "SESTD1"   "HNRNPAB"  "SCMH1"    "JUNB"     "SYTL2"   
## [19] "AGPAT4"   "C12orf75" "TENT5C"   "PKM"      "CD38"     "DISC1"   
## [25] "ST8SIA1" 
## 
## $turquoise
##  [1] "TUBA1B"   "TUBB"     "RRM2"     "STMN1"    "DIAPH3"   "POLQ"    
##  [7] "NCAPG"    "HIST1H1B" "NUSAP1"   "PCLAF"    "ASPM"     "KIF11"   
## [13] "NCAPG2"   "TOP2A"    "CIT"      "KIF15"    "ATAD5"    "KNL1"    
## [19] "BRIP1"    "MKI67"    "CIP2A"    "FANCI"    "STIL"     "ECT2"    
## [25] "SSBP2"    "MCF2L2"   "CD226"    "CENPF"    "CENPE"    "AFF3"    
## [31] "ANTXR2"   "LTB"      "APOLD1"   "DTL"      "FANCA"    "CLSPN"   
## [37] "MCM4"     "HELLS"    "CENPP"    "WDR76"    "FHIT"     "ANXA2"   
## [43] "TGFBR3"   "MCM5"     "MCM7"     "POLA1"    "BRCA2"    "KNTC1"   
## 
## $blue
##  [1] "RPL28"  "RPL18A" "RPL41"  "RPL13A" "RPS15"  "RPL27A" "RPS19"  "RPL11" 
##  [9] "RPS12"  "RPS28"  "RPS23"  "RPS29"  "RPS18"  "RPL10"  "RPL37"  "RPS2"  
## [17] "RPLP1"  "B2M"    "RPS16"  "RPL18"  "RPL37A" "RPS27A" "RPLP2"  "RPS15A"
## [25] "RPL3"  
## 
## [1] "7604/10000 permutations failed the Mann-Whitney test."
## [1] "0/10000 permutations failed the Mann-Whitney test."
## [1] "0/10000 permutations failed the Mann-Whitney test."
## $turquoise
##  [1] "TUBA1B"   "TUBB"     "RRM2"     "STMN1"    "DIAPH3"   "POLQ"    
##  [7] "NCAPG"    "HIST1H1B" "NUSAP1"   "PCLAF"    "ASPM"     "KIF11"   
## [13] "NCAPG2"   "TOP2A"    "CIT"      "KIF15"    "ATAD5"    "KNL1"    
## [19] "BRIP1"    "MKI67"    "CIP2A"    "FANCI"    "STIL"     "ECT2"    
## [25] "SSBP2"    "MCF2L2"   "CD226"    "CENPF"    "CENPE"    "AFF3"    
## [31] "ANTXR2"   "LTB"      "APOLD1"   "DTL"      "FANCA"    "CLSPN"   
## [37] "MCM4"     "HELLS"    "CENPP"    "WDR76"    "FHIT"     "ANXA2"   
## [43] "TGFBR3"   "MCM5"     "MCM7"     "POLA1"    "BRCA2"    "KNTC1"   
## 
## $blue
##  [1] "RPL28"  "RPL18A" "RPL41"  "RPL13A" "RPS15"  "RPL27A" "RPS19"  "RPL11" 
##  [9] "RPS12"  "RPS28"  "RPS23"  "RPS29"  "RPS18"  "RPL10"  "RPL37"  "RPS2"  
## [17] "RPLP1"  "B2M"    "RPS16"  "RPL18"  "RPL37A" "RPS27A" "RPLP2"  "RPS15A"
## [25] "RPL3"

ID modules without merging similar

Rational here is that we might be losing interesting variation by merging. Reality was that we ended up with more modules but they mapped back to the same things.

ModID<-function(seuratobj, cluster=NA, prefix="modules", nPCS="", softpower="", combat=FALSE, batch=NA){
    #subset to a cluster or clusters of interest if desired 
    if (sum(!is.na(cluster))>1){test_clus<-subset(seuratobj, idents=cluster)}else{test_clus<-seuratobj}
    
    #perform feature selection, scaling, and PCA 
    if(combat==TRUE){
        test_clus_2<-test_clus@assays$RNA@data
        test_clus_2<-ComBat(test_clus_2,batch)
        test_clus@assays$RNA@data <- test_clus_2
    }
    test_clus<-FindVariableFeatures(test_clus, verbose=FALSE)
    test_clus<-ScaleData(test_clus, verbose=FALSE)
    test_clus<-RunPCA(test_clus, verbose=FALSE)
    
    #prompot for number of PCs needed 
    if(nPCS==""){print(ElbowPlot(test_clus))}
    while(nPCS==""){nPCS<-readline("how many NPCs?")}
    nPCS<-as.integer(nPCS)
    #expanding the genes by projecting the dims to more PCs 
    test_clus<-ProjectDim(test_clus)
    
    #gere we'll grab the top 25 features (positive and negative) associated with each dimension 
    genes<-c()
    for (i in 1:nPCS){genelist<-TopFeatures(test_clus, dim = i, nfeatures = 50,balanced = T, projected = T )
    for (i in 1:length(genelist)){genes<-c(genes, genelist[[i]])}}
    #remove duplicates 
    genes<-unique(genes)
    #grab the matrix from the seurat object 
    test_clus<-as.matrix(test_clus@assays$RNA@data[genes,])
    
    #prompt for softpower 
    FindPower(datExpr=t(test_clus))
    while(softpower==""){softpower<-readline("what softpower?")}
    softpower<-as.integer(softpower)
    #run the remaining functions to actually cut the tree
    test_clus_tom<-ClusterTOM(datExpr = t(test_clus), softPower = softpower)
    test_clus_mods<-CutTOMTree(datExpr = t(test_clus), geneTree = test_clus_tom$geneTree, dissTOM = test_clus_tom$dissTOM, minModuleSize = 10 )
    #omiting those that will merge modules so that we can test more individual modules to possibly discover more nuanced differences
    test_clus_merge_mods.isSig = sapply(test_clus_mods$modules, function(module){
        TestModuleSignificance(mod = module, dissTOM = test_clus_tom$dissTOM, expr.data = test_clus,
                               n_perm = 10000, pval = 0.05, n.bin = 10)
    })
    test_clus_merge_mods.isSig = test_clus_mods$modules[test_clus_merge_mods.isSig]
    print(test_clus_merge_mods.isSig)
    names(test_clus_merge_mods.isSig)<-paste(prefix,names(test_clus_merge_mods.isSig ), sep="_")
    return(test_clus_merge_mods.isSig)
}

#tried npcs/softpower 

#7/18 - 3 modules identified - proliferation, ribosomal, and second proliferation 
#8/10 - 5 modules, ribosomal + 4 proliferating 
#5/5 - 3 modules identified - proliferation, ribosomal, and second proliferation 

mods<-ModID(results, combat = TRUE, batch = results$batch, nPCS = 7, softpower = 18)
## Found 24380 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.
## Found3batches
## Adjusting for0covariate(s) or covariate level(s)
## Standardizing Data across genes
## Fitting L/S model and finding priors
## Finding parametric adjustments
## Adjusting the Data
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at -2.4183
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.30103
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 9.8548e-16
## PC_ 1 
## Positive:  ANK3, ZBTB20, PCAT1, FAAH2, TSHZ2, MAML2, ABLIM1, PRKCA, SESN3, TTN 
##     LDLRAD4, UST, MPP7, FMN1, AL353660.1, ACYP2, TAFA1, MGAT4A, AL136456.1, EPHA4 
## Negative:  TUBA1B, TUBB, HIST1H4C, RRM2, STMN1, DIAPH3, HIST1H3B, GAPDH, POLQ, ACTB 
##     HIST1H2AH, NCAPG, TK1, HIST1H1B, PFN1, TMPO, NUSAP1, PCLAF, SHCBP1, ASPM 
## PC_ 2 
## Positive:  POLQ, NCAPG, SPC25, DIAPH3, CIT, KIF15, KIF18B, RRM2, ASPM, KIF11 
##     CDCA2, NCAPG2, ATAD5, NUSAP1, KIF14, KIF4A, TOP2A, BUB1B, KNL1, BRIP1 
## Negative:  RPL28, RPL18A, RPL41, RPL13A, RPS15, RPL27A, RPS19, RPL11, RPS12, RPS28 
##     RPS23, RPS29, RPS18, RPL10, RPL37, RPS27, RPS2, RPLP1, B2M, RPS16 
## PC_ 3 
## Positive:  ABLIM1, MGAT5, RORA, MAP3K1, ZBTB20, SSBP2, WDFY2, MAP3K4, MCF2L2, ATP13A3 
##     AL136456.1, RFX3, ARL15, TSHZ2, PHF21A, STAM, DUSP16, AKT3, MAST4, UST 
## Negative:  MTRNR2L8, NKG7, AGAP1, CCL5, PECAM1, CD8A, MKI67, CENPF, GZMA, KLRK1 
##     HMGB2, TOP2A, CD8B, INCENP, KPNA2, AOAH, HIST1H1E, CENPE, CDCA2, DLGAP5 
## PC_ 4 
## Positive:  PECAM1, PTPN3, DTHD1, MIR4422HG, LYST, MCTP2, NCOA1, GZMK, KLRK1, CNR2 
##     SESTD1, NKG7, AGAP1, AL031599.1, HAVCR2, KCNQ5, UBE2E2, SAMD3, GZMA, AFF3 
## Negative:  LTB, MAP3K1, SORL1, KIF14, KIF4A, CDCA2, CDCA8, UBE2C, KIF23, TSHZ2 
##     ANK3, CKAP2L, KIF2C, PLK1, KIF20A, CCNB2, ESPL1, CCNF, GTSE1, HMMR 
## PC_ 5 
## Positive:  MIR4422HG, PTPN3, SESTD1, DTHD1, AFF3, GZMK, SESN3, AL031599.1, LDLRAD4, PECAM1 
##     UBE2E2, ANTXR2, BACH2, LYST, HAVCR2, KIF20A, PLK1, DSCAML1, PACSIN1, KIF14 
## Negative:  EXO1, DTL, FANCA, CLSPN, TRERF1, CCNE2, CDC6, MSH6, CDC45, CTSW 
##     MCM4, POLE2, FAM111B, HELLS, MCOLN2, CENPP, C1orf21, CCL5, MPP6, MCM10 
##    Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
## 1      1  0.00209  0.443         0.1190 136.000  1.34e+02 151.00
## 2      2  0.23500 -2.250         0.3650  73.400  7.07e+01  94.00
## 3      3  0.53100 -2.360         0.5380  40.800  3.77e+01  61.60
## 4      4  0.64600 -2.090         0.5780  23.400  2.03e+01  42.40
## 5      5  0.66600 -1.850         0.5830  13.900  1.11e+01  30.40
## 6      6  0.60200 -1.900         0.5280   8.620  6.09e+00  22.60
## 7      7  0.60300 -1.620         0.5240   5.540  3.41e+00  17.20
## 8      8  0.62800 -1.480         0.5540   3.700  1.93e+00  13.40
## 9      9  0.66600 -1.390         0.6050   2.550  1.14e+00  10.60
## 10    10  0.67600 -1.310         0.6300   1.820  6.62e-01   8.50
## 11    12  0.21200 -2.540         0.1180   0.989  2.33e-01   5.59
## 12    14  0.20500 -2.180         0.1330   0.578  9.50e-02   3.79
## 13    16  0.21000 -2.100         0.1550   0.354  4.05e-02   2.64
## 14    18  0.80100 -1.090         0.8340   0.223  1.68e-02   1.86
## 15    20  0.19200 -1.870         0.0244   0.144  7.77e-03   1.33

## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.

## dynamicMods
##   0   1   2   3   4 
## 165  33  25  22  13

## [1] "4129/10000 permutations failed the Mann-Whitney test."
## [1] "0/10000 permutations failed the Mann-Whitney test."
## [1] "206/10000 permutations failed the Mann-Whitney test."
## [1] "0/10000 permutations failed the Mann-Whitney test."
## $turquoise
##  [1] "TUBA1B"   "TUBB"     "STMN1"    "DIAPH3"   "POLQ"     "NCAPG"   
##  [7] "HIST1H1B" "PCLAF"    "ASPM"     "KIF11"    "NCAPG2"   "TOP2A"   
## [13] "CIT"      "KIF15"    "ATAD5"    "KNL1"     "CIP2A"    "STIL"    
## [19] "SSBP2"    "CD226"    "CENPF"    "ANTXR2"   "LTB"      "DTL"     
## [25] "FANCA"    "HELLS"    "WDR76"    "FHIT"     "ANXA2"    "TGFBR3"  
## [31] "MCM5"     "MCM7"     "POLA1"   
## 
## $yellow
##  [1] "RRM2"   "NUSAP1" "BRIP1"  "MKI67"  "FANCI"  "ECT2"   "MCF2L2" "CENPE" 
##  [9] "AFF3"   "APOLD1" "CLSPN"  "MCM4"   "CENPP" 
## 
## $blue
##  [1] "RPL28"  "RPL18A" "RPL41"  "RPL13A" "RPS15"  "RPL27A" "RPS19"  "RPL11" 
##  [9] "RPS12"  "RPS28"  "RPS23"  "RPS29"  "RPS18"  "RPL10"  "RPL37"  "RPS27" 
## [17] "RPS2"   "RPLP1"  "B2M"    "RPS16"  "RPL18"  "RPL37A" "RPS27A" "RPLP2" 
## [25] "RPS15A"
devtools::session_info()
## Warning in system("timedatectl", intern = TRUE): running command 'timedatectl'
## had status 1
## - Session info ---------------------------------------------------------------
##  setting  value
##  version  R version 4.2.0 (2022-04-22)
##  os       Red Hat Enterprise Linux 8.7 (Ootpa)
##  system   x86_64, linux-gnu
##  ui       X11
##  language (EN)
##  collate  C
##  ctype    C
##  tz       Etc/UTC
##  date     2023-06-01
##  pandoc   2.17.1.1 @ /usr/lib/rstudio-server/bin/quarto/bin/ (via rmarkdown)
## 
## - Packages -------------------------------------------------------------------
##  package          * version   date (UTC) lib source
##  abind              1.4-5     2016-07-21 [2] CRAN (R 4.2.0)
##  annotate           1.76.0    2022-11-01 [1] Bioconductor
##  AnnotationDbi      1.60.2    2023-03-10 [1] Bioconductor
##  backports          1.4.1     2021-12-13 [2] CRAN (R 4.2.0)
##  base64enc          0.1-3     2015-07-28 [2] CRAN (R 4.2.0)
##  beeswarm           0.4.0     2021-06-01 [2] CRAN (R 4.2.0)
##  Biobase            2.58.0    2022-11-01 [1] Bioconductor
##  BiocGenerics       0.44.0    2022-11-01 [1] Bioconductor
##  BiocParallel     * 1.32.6    2023-03-17 [1] Bioconductor
##  Biostrings         2.66.0    2022-11-01 [1] Bioconductor
##  bit                4.0.5     2022-11-15 [2] CRAN (R 4.2.0)
##  bit64              4.0.5     2020-08-30 [2] CRAN (R 4.2.0)
##  bitops             1.0-7     2021-04-24 [2] CRAN (R 4.2.0)
##  blob               1.2.4     2023-03-17 [1] CRAN (R 4.2.0)
##  bslib              0.4.2     2022-12-16 [1] CRAN (R 4.2.0)
##  cachem             1.0.8     2023-05-01 [1] CRAN (R 4.2.0)
##  callr              3.7.3     2022-11-02 [1] CRAN (R 4.2.0)
##  checkmate          2.1.0     2022-04-21 [2] CRAN (R 4.2.0)
##  cli                3.6.1     2023-03-23 [1] CRAN (R 4.2.0)
##  cluster            2.1.4     2022-08-22 [2] CRAN (R 4.2.0)
##  codetools          0.2-19    2023-02-01 [2] CRAN (R 4.2.0)
##  colorspace         2.1-0     2023-01-23 [2] CRAN (R 4.2.0)
##  cowplot          * 1.1.1     2020-12-30 [2] CRAN (R 4.2.0)
##  crayon             1.5.2     2022-09-29 [2] CRAN (R 4.2.0)
##  data.table         1.14.8    2023-02-17 [2] CRAN (R 4.2.0)
##  DBI                1.1.3     2022-06-18 [2] CRAN (R 4.2.0)
##  deldir             1.0-6     2021-10-23 [2] CRAN (R 4.2.0)
##  devtools           2.4.5     2022-10-11 [1] CRAN (R 4.2.0)
##  digest             0.6.31    2022-12-11 [2] CRAN (R 4.2.0)
##  doParallel       * 1.0.17    2022-02-07 [1] CRAN (R 4.2.0)
##  dplyr            * 1.1.2     2023-04-20 [1] CRAN (R 4.2.0)
##  dynamicTreeCut   * 1.63-1    2016-03-11 [2] CRAN (R 4.2.0)
##  edgeR              3.40.2    2023-01-19 [1] Bioconductor
##  ellipsis           0.3.2     2021-04-29 [2] CRAN (R 4.2.0)
##  evaluate           0.20      2023-01-17 [2] CRAN (R 4.2.0)
##  fansi              1.0.4     2023-01-22 [2] CRAN (R 4.2.0)
##  farver             2.1.1     2022-07-06 [2] CRAN (R 4.2.0)
##  fastcluster      * 1.2.3     2021-05-24 [2] CRAN (R 4.2.0)
##  fastmap            1.1.1     2023-02-24 [1] CRAN (R 4.2.0)
##  fastmatch          1.1-3     2021-07-23 [2] CRAN (R 4.2.0)
##  fitdistrplus       1.1-8     2022-03-10 [2] CRAN (R 4.2.0)
##  flashClust       * 1.01-2    2012-08-21 [1] CRAN (R 4.2.0)
##  foreach          * 1.5.2     2022-02-02 [2] CRAN (R 4.2.0)
##  foreign            0.8-84    2022-12-06 [2] CRAN (R 4.2.0)
##  Formula            1.2-5     2023-02-24 [1] CRAN (R 4.2.0)
##  fs                 1.6.1     2023-02-06 [2] CRAN (R 4.2.0)
##  future             1.32.0    2023-03-07 [1] CRAN (R 4.2.0)
##  future.apply       1.10.0    2022-11-05 [1] CRAN (R 4.2.0)
##  genefilter       * 1.80.3    2023-01-19 [1] Bioconductor
##  generics           0.1.3     2022-07-05 [2] CRAN (R 4.2.0)
##  GenomeInfoDb       1.34.9    2023-02-02 [1] Bioconductor
##  GenomeInfoDbData   1.2.9     2023-03-17 [1] Bioconductor
##  GenomicRanges      1.50.2    2022-12-16 [1] Bioconductor
##  ggbeeswarm         0.7.2     2023-04-29 [1] CRAN (R 4.2.0)
##  ggplot2          * 3.4.2     2023-04-03 [1] CRAN (R 4.2.0)
##  ggrastr            1.0.1     2021-12-08 [1] CRAN (R 4.2.0)
##  ggrepel            0.9.3     2023-02-03 [1] CRAN (R 4.2.0)
##  ggridges           0.5.4     2022-09-26 [1] CRAN (R 4.2.0)
##  globals            0.16.2    2022-11-21 [1] CRAN (R 4.2.0)
##  glue               1.6.2     2022-02-24 [2] CRAN (R 4.2.0)
##  GO.db              3.16.0    2023-03-17 [1] Bioconductor
##  goftest            1.2-3     2021-10-07 [2] CRAN (R 4.2.0)
##  gridExtra          2.3       2017-09-09 [2] CRAN (R 4.2.0)
##  gtable             0.3.3     2023-03-21 [1] CRAN (R 4.2.0)
##  highr              0.10      2022-12-22 [1] CRAN (R 4.2.0)
##  Hmisc            * 5.1-0     2023-05-08 [1] CRAN (R 4.2.0)
##  htmlTable          2.4.1     2022-07-07 [1] CRAN (R 4.2.0)
##  htmltools          0.5.5     2023-03-23 [1] CRAN (R 4.2.0)
##  htmlwidgets        1.6.2     2023-03-17 [1] CRAN (R 4.2.0)
##  httpuv             1.6.9     2023-02-14 [1] CRAN (R 4.2.0)
##  httr               1.4.5     2023-02-24 [1] CRAN (R 4.2.0)
##  ica                1.0-3     2022-07-08 [2] CRAN (R 4.2.0)
##  igraph             1.4.2     2023-04-07 [1] CRAN (R 4.2.0)
##  impute             1.72.3    2023-01-19 [1] Bioconductor
##  IRanges            2.32.0    2022-11-01 [1] Bioconductor
##  irlba              2.3.5.1   2022-10-03 [1] CRAN (R 4.2.0)
##  iterators        * 1.0.14    2022-02-05 [2] CRAN (R 4.2.0)
##  jquerylib          0.1.4     2021-04-26 [2] CRAN (R 4.2.0)
##  jsonlite           1.8.4     2022-12-06 [2] CRAN (R 4.2.0)
##  KEGGREST           1.38.0    2022-11-01 [1] Bioconductor
##  KernSmooth         2.23-20   2021-05-03 [2] CRAN (R 4.2.0)
##  knitr              1.42      2023-01-25 [1] CRAN (R 4.2.0)
##  labeling           0.4.2     2020-10-20 [2] CRAN (R 4.2.0)
##  later              1.3.0     2021-08-18 [2] CRAN (R 4.2.0)
##  lattice            0.21-8    2023-04-05 [1] CRAN (R 4.2.0)
##  lazyeval           0.2.2     2019-03-15 [2] CRAN (R 4.2.0)
##  leiden             0.4.3     2022-09-10 [1] CRAN (R 4.2.0)
##  lifecycle          1.0.3     2022-10-07 [1] CRAN (R 4.2.0)
##  limma              3.54.2    2023-02-28 [1] Bioconductor
##  listenv            0.9.0     2022-12-16 [2] CRAN (R 4.2.0)
##  lmtest             0.9-40    2022-03-21 [2] CRAN (R 4.2.0)
##  locfit             1.5-9.7   2023-01-02 [2] CRAN (R 4.2.0)
##  magrittr           2.0.3     2022-03-30 [2] CRAN (R 4.2.0)
##  MASS               7.3-59    2023-04-21 [1] CRAN (R 4.2.0)
##  Matrix             1.5-4     2023-04-04 [1] CRAN (R 4.2.0)
##  matrixStats        0.63.0    2022-11-18 [2] CRAN (R 4.2.0)
##  memoise            2.0.1     2021-11-26 [2] CRAN (R 4.2.0)
##  mgcv             * 1.8-42    2023-03-02 [1] CRAN (R 4.2.0)
##  mime               0.12      2021-09-28 [2] CRAN (R 4.2.0)
##  miniUI             0.1.1.1   2018-05-18 [2] CRAN (R 4.2.0)
##  munsell            0.5.0     2018-06-12 [2] CRAN (R 4.2.0)
##  nlme             * 3.1-162   2023-01-31 [1] CRAN (R 4.2.0)
##  nnet               7.3-19    2023-05-03 [1] CRAN (R 4.2.0)
##  openxlsx         * 4.2.5.2   2023-02-06 [1] CRAN (R 4.2.0)
##  parallelly         1.35.0    2023-03-23 [1] CRAN (R 4.2.0)
##  patchwork          1.1.2     2022-08-19 [1] CRAN (R 4.2.0)
##  pbapply            1.7-0     2023-01-13 [1] CRAN (R 4.2.0)
##  pillar             1.9.0     2023-03-22 [1] CRAN (R 4.2.0)
##  pkgbuild           1.4.0     2022-11-27 [1] CRAN (R 4.2.0)
##  pkgconfig          2.0.3     2019-09-22 [2] CRAN (R 4.2.0)
##  pkgload            1.3.2     2022-11-16 [1] CRAN (R 4.2.0)
##  plotly             4.10.1    2022-11-07 [1] CRAN (R 4.2.0)
##  plyr               1.8.8     2022-11-11 [1] CRAN (R 4.2.0)
##  png                0.1-8     2022-11-29 [1] CRAN (R 4.2.0)
##  polyclip           1.10-4    2022-10-20 [1] CRAN (R 4.2.0)
##  preprocessCore     1.60.2    2023-01-19 [1] Bioconductor
##  prettyunits        1.1.1     2020-01-24 [2] CRAN (R 4.2.0)
##  processx           3.8.1     2023-04-18 [1] CRAN (R 4.2.0)
##  profvis            0.3.8     2023-05-02 [1] CRAN (R 4.2.0)
##  progressr          0.13.0    2023-01-10 [1] CRAN (R 4.2.0)
##  promises           1.2.0.1   2021-02-11 [2] CRAN (R 4.2.0)
##  ps                 1.7.5     2023-04-18 [1] CRAN (R 4.2.0)
##  purrr              1.0.1     2023-01-10 [1] CRAN (R 4.2.0)
##  R6                 2.5.1     2021-08-19 [2] CRAN (R 4.2.0)
##  RANN               2.6.1     2019-01-08 [2] CRAN (R 4.2.0)
##  RColorBrewer       1.1-3     2022-04-03 [2] CRAN (R 4.2.0)
##  Rcpp               1.0.10    2023-01-22 [1] CRAN (R 4.2.0)
##  RcppAnnoy          0.0.20    2022-10-27 [1] CRAN (R 4.2.0)
##  RcppRoll           0.3.0     2018-06-05 [2] CRAN (R 4.2.0)
##  RCurl              1.98-1.12 2023-03-27 [1] CRAN (R 4.2.0)
##  remotes            2.4.2     2021-11-30 [2] CRAN (R 4.2.0)
##  reshape2           1.4.4     2020-04-09 [2] CRAN (R 4.2.0)
##  reticulate         1.28      2023-01-27 [1] CRAN (R 4.2.0)
##  rlang              1.1.1     2023-04-28 [1] CRAN (R 4.2.0)
##  rmarkdown          2.21      2023-03-26 [1] CRAN (R 4.2.0)
##  ROCR               1.0-11    2020-05-02 [2] CRAN (R 4.2.0)
##  rpart              4.1.19    2022-10-21 [1] CRAN (R 4.2.0)
##  Rsamtools          2.14.0    2022-11-01 [1] Bioconductor
##  RSQLite            2.3.1     2023-04-03 [1] CRAN (R 4.2.0)
##  rstudioapi         0.14      2022-08-22 [1] CRAN (R 4.2.0)
##  Rtsne              0.16      2022-04-17 [2] CRAN (R 4.2.0)
##  S4Vectors          0.36.2    2023-02-26 [1] Bioconductor
##  sass               0.4.5     2023-01-24 [1] CRAN (R 4.2.0)
##  scales             1.2.1     2022-08-20 [1] CRAN (R 4.2.0)
##  scattermore        0.8       2022-02-14 [1] CRAN (R 4.2.0)
##  sctransform        0.3.5     2022-09-21 [1] CRAN (R 4.2.0)
##  sessioninfo        1.2.2     2021-12-06 [2] CRAN (R 4.2.0)
##  Seurat           * 4.3.0     2022-11-18 [1] CRAN (R 4.2.0)
##  SeuratObject     * 4.1.3     2022-11-07 [1] CRAN (R 4.2.0)
##  shiny              1.7.4     2022-12-15 [1] CRAN (R 4.2.0)
##  Signac           * 1.9.0     2022-12-08 [1] CRAN (R 4.2.0)
##  sp                 1.6-0     2023-01-19 [1] CRAN (R 4.2.0)
##  spatstat.data      3.0-1     2023-03-12 [1] CRAN (R 4.2.0)
##  spatstat.explore   3.1-0     2023-03-14 [1] CRAN (R 4.2.0)
##  spatstat.geom      3.1-0     2023-03-12 [1] CRAN (R 4.2.0)
##  spatstat.random    3.1-4     2023-03-13 [1] CRAN (R 4.2.0)
##  spatstat.sparse    3.0-1     2023-03-12 [1] CRAN (R 4.2.0)
##  spatstat.utils     3.0-2     2023-03-11 [1] CRAN (R 4.2.0)
##  stringi            1.7.12    2023-01-11 [1] CRAN (R 4.2.0)
##  stringr            1.5.0     2022-12-02 [1] CRAN (R 4.2.0)
##  survival           3.5-5     2023-03-12 [1] CRAN (R 4.2.0)
##  sva              * 3.46.0    2022-11-01 [1] Bioconductor
##  tensor             1.5       2012-05-05 [2] CRAN (R 4.2.0)
##  tibble             3.2.1     2023-03-20 [1] CRAN (R 4.2.0)
##  tidyr              1.3.0     2023-01-24 [1] CRAN (R 4.2.0)
##  tidyselect         1.2.0     2022-10-10 [1] CRAN (R 4.2.0)
##  urlchecker         1.0.1     2021-11-30 [1] CRAN (R 4.2.0)
##  usethis            2.1.6     2022-05-25 [1] CRAN (R 4.2.0)
##  utf8               1.2.3     2023-01-31 [1] CRAN (R 4.2.0)
##  uwot               0.1.14    2022-08-22 [1] CRAN (R 4.2.0)
##  vctrs              0.6.2     2023-04-19 [1] CRAN (R 4.2.0)
##  vipor              0.4.5     2017-03-22 [2] CRAN (R 4.2.0)
##  viridisLite        0.4.2     2023-05-02 [1] CRAN (R 4.2.0)
##  WGCNA            * 1.72-1    2023-01-18 [1] CRAN (R 4.2.0)
##  withr              2.5.0     2022-03-03 [2] CRAN (R 4.2.0)
##  xfun               0.39      2023-04-20 [1] CRAN (R 4.2.0)
##  XML                3.99-0.14 2023-03-19 [1] CRAN (R 4.2.0)
##  xtable             1.8-4     2019-04-21 [2] CRAN (R 4.2.0)
##  XVector            0.38.0    2022-11-01 [1] Bioconductor
##  yaml               2.3.7     2023-01-23 [1] CRAN (R 4.2.0)
##  zip                2.3.0     2023-04-17 [1] CRAN (R 4.2.0)
##  zlibbioc           1.44.0    2022-11-01 [1] Bioconductor
##  zoo                1.8-12    2023-04-13 [1] CRAN (R 4.2.0)
## 
##  [1] /gpfs/gibbs/project/ya-chi_ho/jac369/R/4.2
##  [2] /vast/palmer/apps/avx2/software/R/4.2.0-foss-2020b/lib64/R/library
## 
## ------------------------------------------------------------------------------